
##########################################################################

  Table 2. Monte Carlo Simulation Results for Asymptotic Properties

##########################################################################

library(lavaan)

identity<-diag(1, 24,24)
phi<-matrix(c( 1, 0.3, 0.4,
               0.3, 1, 0.5,
               0.4, 0.5, 1
), 3,3)

ld<-matrix(c(0.65, 0.65, 0.7, 0.7, 0.7, 0.6, 0.6, 0.5,   0.6, 0, 0, 0, 0, 0, 0,0,  0, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0.5, 0, 0,  0.6, 0.6, 0.6, 0.7, 0.7, 0.5, 0.5, 0.55, 0, 0, 0, 0.5, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0.55,  0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.7, 0.7 
), nrow=24, ncol=3)

ld2<-t(ld)
psi<-diag(1,24)-diag(ld %*% phi %*% ld2)
diag(psi)




popmodel<-'
f1 =~ 0.65*x1 + 0.65*x2 + 0.7*x3 + 0.7*x4 + 0.7*x5 + 0.7*x6 + 0.6*x7 + 0.5*x8 + 0.5*x9
f2 =~ 0.6*x9 + 0.6*x10 + 0.6*x11 + 0.7*x12 + 0.7*x13 + 0.5*x14 + 0.5*x15 + 0.65*x16 + 0.5*x6 + 0.75*x20
f3 =~ 0.5*x17 + 0.5*x18 + 0.5*x19 + 0.6*x20 + 0.6*x21 + 0.6*x22 + 0.7*x23 + 0.7*x24 + 0.45*x16

f3 ~~ 0.4*f1 + 0.5*f2
f1 ~~ 0.3*f2

f1 ~ 0*1
f2 ~ 0*1
f3 ~ 0*1

x1 ~~ 0.5775*x1
x2 ~~ 0.5775*x2

x3 ~~ 0.51*x3
x4 ~~ 0.51*x4
x5 ~~ 0.51*x5
x6 ~~ 0.05*x6

x7 ~~ 0.64*x7
x8 ~~ 0.21*x8
x9 ~~ 0.21*x9

x10 ~~  0.64*x10
x11 ~~  0.51*x11

x12 ~~ 0.51*x12
x13 ~~ 0.75*x13
x14 ~~ 0.75*x14
x15 ~~ 0.5775*x15

x16 ~~ 0.7975*x16
x17 ~~ 0.75*x17

x18 ~~ 0.75*x18
x19 ~~ 0.25*x19
x20 ~~ 0.64*x20
x21 ~~ 0.64*x21
x22 ~~ 0.64*x22

x23 ~~ 0.51*x23
x24 ~~ 0.51*x24


'

mymodel<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8 + x9
f2 =~ x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16 + x6 + x20
f3 =~ x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x16
f3 ~~ f1
f1 ~~ f2

'


################################

      N=100

###############################
set.seed(100)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=100L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p


sd(bin.1000[, 1], na.rm = TRUE)


df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)


quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

      N=150

###############################
set.seed(150)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=150L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)


p <- mean(bin.1000[, 2], na.rm = TRUE)
p




df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



###############################

      N=200

###############################
set.seed(200)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=200L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)


p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

      N=250

###############################
set.seed(250)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=250L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1


sd(bin.1000[, 1], na.rm = TRUE)


p <- mean(bin.1000[, 2], na.rm = TRUE)
p


df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)


quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)


##############################

      N=300

##############################
set.seed(300)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=300L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

      N=400

##############################
set.seed(400)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=400L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

      N=500

##############################
set.seed(500)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=500L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p


df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)


quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)





##############################

      N=800

#############################

set.seed(800)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=800L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)

quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

      N=1000

##############################


set.seed(1000)

bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=1000L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1


sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p

df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)


quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



##############################

        N=2000

##############################

set.seed(2000)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=2000L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)

quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)



################################

      N=5000

################################

set.seed(5000)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=5000L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)

p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)




###############################

          N=10000

###############################
set.seed(10000)
bin.1000 <- matrix(NA,nrow=1000, ncol=5)
for (i in 1:1000) {
  mydata<-simulateData(popmodel, model.type = "sem", meanstructure = FALSE, sample.nobs=10000L)
  fit1 <- cfa(model = mymodel, data=mydata, meanstructure=FALSE, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  fitted(fit1)[1]
  pvalue <- lavInspect(fit1, "test")[[1]]$pvalue
  stat <- lavInspect(fit1, "test")[[1]]$stat
  nfi <- fitMeasures(fit1, "nfi")
  cfi <- fitMeasures(fit1, "cfi")
  rmsea<- fitMeasures(fit1, "rmsea")
  
  bin.1000[i,1]<-stat
  bin.1000[i,2]<-pvalue
  bin.1000[i,3]<-nfi
  bin.1000[i,4]<-cfi
  bin.1000[i,5]<-rmsea
  
}


T1 <- mean(bin.1000[, 1], na.rm = TRUE)
T1

sd(bin.1000[, 1], na.rm = TRUE)


p <- mean(bin.1000[, 2], na.rm = TRUE)
p



df<- lavInspect(fit1, "test")[[1]]$df


library(plyr)
x<-sum(bin.1000[,2]<0.05)
total<-length(bin.1000[,2])
rej_rate<-x/total
print(rej_rate)



quantile(bin.1000[, 3],.025)
quantile(bin.1000 [, 3],.975)

quantile(bin.1000[, 4],.025)
quantile(bin.1000 [, 4],.975)

quantile(bin.1000[, 5],.025)
quantile(bin.1000 [, 5],.975)




#######. The End

###### All results match the manuscript, 12/5/24
